library(data.table)
library(magrittr)
library(parallel)
library(dplyr)
library(boot)
library(INLA)

setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM")
dir1="./bayes_pvalue_beta0/"
dir2="./bayes_BF/"
file=read.csv("./20201120/at.least.one.AShM.in.DC.add.BF.beta0.csv",head=T)
file1=file[file$pattern.not.rm.dupl.num.DC>1,]
id=as.character(file$unitID)
group1=c("X2B_X1T","M8_M7","M6_M5","M2_M1","M48_M47","M50_M49","M28_M27","M30_M29","M26_M25","M35_M36","M18_M17","M20_M19","M22_M21","M40_M39")
i=7
fn1=paste0(dir1,group1[i],".bayes_p.txt")
file1=read.table(fn1,head=T,sep = "\t")
file1$unitID=paste(file1$chrom,file1$position,file1$ref,file1$var,sep=":")
file1=data.frame(unitID=file1$unitID,normal_reads1=file1$normal_reads1,normal_reads2=file1$normal_reads2,tumor_reads1=file1$tumor_reads1,tumor_reads2=file1$tumor_reads2)
file1=file1[file1$unitID %in% id,]
names(file1)=c("unitID",paste(group1[i],names(file1)[2:5],sep="_"))
file=file1
for(i in 8:10){
  fn1=paste0(dir1,group1[i],".bayes_p.txt")
  file1=read.table(fn1,head=T,sep = "\t")
  file1$unitID=paste(file1$chrom,file1$position,file1$ref,file1$var,sep=":")
  file1=data.frame(unitID=file1$unitID,normal_reads1=file1$normal_reads1,normal_reads2=file1$normal_reads2,tumor_reads1=file1$tumor_reads1,tumor_reads2=file1$tumor_reads2)
  file1=file1[file1$unitID %in% id,]
  names(file1)=c("unitID",paste(group1[i],names(file1)[2:5],sep="_"))
  file=merge(file,file1,by="unitID",all=T)
}

filekp=file
refnum=grep(names(filekp),pattern="reads1")###左边是ref，右边是alt
altnum=grep(names(filekp),pattern="reads2")
file=filekp[,c(refnum,altnum)]	###control
result=matrix(,nrow(file),ncol=4)
for (i in 1:nrow(file)){
  reads=as.numeric(as.matrix(file)[i,])
  ref=reads[1:(ncol(file)/2)]
  var=reads[(ncol(file)/2+1):ncol(file)]
  tid=rep((1:4)-1,each=2)
  sl=which(!is.na(ref))
  if ((length(sl)/2)>=1){
    df=data.frame(y=var[sl],Totalreads=ref[sl]+var[sl],x2=tid[sl])
    null <- logit(median(1-df[,1]/df[,2]))
    formula = y ~ 1 + f(x2, model = "iid")
    #formula = ALT_COUNT ~ 1 + f(TISSUE_ID, model = "iid")
    m1 <- inla(formula, data = df, family = "binomial", Ntrials = Totalreads,quantile = c(0.005, 0.025, 0.975, 0.995))
    m <- m1$marginals.fixed[[1]]
    lower_p <- inla.pmarginal(null, m)
    upper_p <- 1 - inla.pmarginal(null, m)
    post_pred_p <- 2 * (min(lower_p, upper_p))
    coef <- m1$summary.fixed
    result[i,1]=post_pred_p
    result[i,2]=coef$mean
    result[i,3]=coef$`0.025quant`
    result[i,4]=coef$`0.975quant`
  }else{
    df=data.frame(y=var[sl],Totalreads=ref[sl]+var[sl])
    null <- logit(median(1-df[,1]/df[,2]))
    formula = y ~ 1
    m1 <- inla(formula, data = df, family = "binomial", Ntrials = Totalreads,quantile = c(0.005, 0.025, 0.975, 0.995))
    m <- m1$marginals.fixed[[1]]
    lower_p <- inla.pmarginal(null, m)
    upper_p <- 1 - inla.pmarginal(null, m)
    post_pred_p <- 2 * (min(lower_p, upper_p))
    coef <- m1$summary.fixed
    #ci95 <- c(coef[4], coef[5])
    result[i,1]=post_pred_p
    result[i,2]=coef$mean
    result[i,3]=coef$`0.025quant`
    result[i,4]=coef$`0.975quant`
  }
}
filekp$CC.pvalue=result[,1]
filekp$CC.FDR=p.adjust(filekp$CC.pvalue,method = "BH")
filekp$CC.beta0=result[,2]
filekp$CC.beta0.lower=result[,3]
filekp$CC.beta0.upper=result[,4]
rt=data.frame(unitID=filekp$unitID,filekp[,18:22])
write.table(rt,"./20201120/CC.beta0.txt",quote=F,row.names = F,sep = "\t")

###健康的
i=11
fn1=paste0(dir1,group1[i],".bayes_p.txt")
file1=read.table(fn1,head=T,sep = "\t")
file1$unitID=paste(file1$chrom,file1$position,file1$ref,file1$var,sep=":")
file1=data.frame(unitID=file1$unitID,normal_reads1=file1$normal_reads1,normal_reads2=file1$normal_reads2,tumor_reads1=file1$tumor_reads1,tumor_reads2=file1$tumor_reads2)
file1=file1[file1$unitID %in% id,]
names(file1)=c("unitID",paste(group1[i],names(file1)[2:5],sep="_"))
file=file1
for(i in 12:14){
  fn1=paste0(dir1,group1[i],".bayes_p.txt")
  file1=read.table(fn1,head=T,sep = "\t")
  file1$unitID=paste(file1$chrom,file1$position,file1$ref,file1$var,sep=":")
  file1=data.frame(unitID=file1$unitID,normal_reads1=file1$normal_reads1,normal_reads2=file1$normal_reads2,tumor_reads1=file1$tumor_reads1,tumor_reads2=file1$tumor_reads2)
  file1=file1[file1$unitID %in% id,]
  names(file1)=c("unitID",paste(group1[i],names(file1)[2:5],sep="_"))
  file=merge(file,file1,by="unitID",all=T)
}

filekp=file
refnum=grep(names(filekp),pattern="reads1")###左边是ref，右边是alt
altnum=grep(names(filekp),pattern="reads2")
file=filekp[,c(refnum,altnum)]	###control
result=matrix(,nrow(file),ncol=4)
for (i in 1:nrow(file)){
  reads=as.numeric(as.matrix(file)[i,])
  ref=reads[1:(ncol(file)/2)]
  var=reads[(ncol(file)/2+1):ncol(file)]
  tid=rep((1:4)-1,each=2)
  sl=which(!is.na(ref))
  if ((length(sl)/2)>=1){
    df=data.frame(y=var[sl],Totalreads=ref[sl]+var[sl],x2=tid[sl])
    null <- logit(median(1-df[,1]/df[,2]))
    formula = y ~ 1 + f(x2, model = "iid")
    #formula = ALT_COUNT ~ 1 + f(TISSUE_ID, model = "iid")
    m1 <- inla(formula, data = df, family = "binomial", Ntrials = Totalreads,quantile = c(0.005, 0.025, 0.975, 0.995))
    m <- m1$marginals.fixed[[1]]
    lower_p <- inla.pmarginal(null, m)
    upper_p <- 1 - inla.pmarginal(null, m)
    post_pred_p <- 2 * (min(lower_p, upper_p))
    coef <- m1$summary.fixed
    result[i,1]=post_pred_p
    result[i,2]=coef$mean
    result[i,3]=coef$`0.025quant`
    result[i,4]=coef$`0.975quant`
  }else{
    df=data.frame(y=var[sl],Totalreads=ref[sl]+var[sl])
    null <- logit(median(1-df[,1]/df[,2]))
    formula = y ~ 1
    m1 <- inla(formula, data = df, family = "binomial", Ntrials = Totalreads,quantile = c(0.005, 0.025, 0.975, 0.995))
    m <- m1$marginals.fixed[[1]]
    lower_p <- inla.pmarginal(null, m)
    upper_p <- 1 - inla.pmarginal(null, m)
    post_pred_p <- 2 * (min(lower_p, upper_p))
    coef <- m1$summary.fixed
    #ci95 <- c(coef[4], coef[5])
    result[i,1]=post_pred_p
    result[i,2]=coef$mean
    result[i,3]=coef$`0.025quant`
    result[i,4]=coef$`0.975quant`
  }
}
filekp$HC.pvalue=result[,1]
filekp$HC.FDR=p.adjust(filekp$HC.pvalue,method = "BH")
filekp$HC.beta0=result[,2]
filekp$HC.beta0.lower=result[,3]
filekp$HC.beta0.upper=result[,4]
rt=data.frame(unitID=filekp$unitID,filekp[,18:22])
write.table(rt,"./20201120/HC.beta0.txt",quote=F,row.names = F,sep = "\t")